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ABSTRACT 


The estimation of time delay and Doppler difference of a signal arriving at two 
physically separated sensors is investigated in this thesis. Usually, modified cross power 
spectrum coupled with Doppler compensation is used to detect a common, passive signal 
received at two separated sensors. Another successful approach uses the cross coherence 
to achieve this goal. This thesis modifies these two techniques to model the Doppler 
difference via an autoregressive (AR) technique. Analytical results are derived and ex- 
perimentallv verified via a computer simulation. Performance at high and low signal to 


noise ratios { S.\ As ) is examined. 
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I. INTRODUCTION 


This thesis investigates the use of autoregressive (AR) models for estimating the 
Doppler difference and differential ttme delay by processing of a narrow band signal 
enutted from a moving source and received at two physically separated sensors. If the 
signal is received at two different geographical positions even in the presence of uncor- 
related noise. then. depending on the signal strength and duration, it 1s possible to esti- 
mate the differential time delav. 

Because the source is moving, the signals that are received at the sensors mav have 
different frequencies due to the Doppler effect. To obtain accurate differential time de- 
lav estimation. Doppler difference compensation is usually required. 

This compensation can be implemented by using frequency shifting of the narrow 
band components of the received signal. This frequency shift can be obtained using a 
Fourier transform. In this thesis we use an AR model to detect the frequency shift. 
Using this Doppler compensation, an estimate of differential time delay can be obtained. 
Estimating the delav and Doppler using an AR model can be interpreted as a form of a 
narrow band coherence procedure, provided the estimations are properly normalized. 

This thesis 1s arranged in six chapters and six appendices. Because the estimation 
of the time delay and Doppler is intimately related to the coherence between two trans- 
formed complex signals. an extensive investigation of coherence is given in Chapter II. 
In Chapter II}. AR models. advantages of AR modeling. and AR model order selection 
are presented. Chapter IV is devoted to the analvsis and the processing of noisy signals 
to estimate the differential time delay and the Doppler difference. To estimate the dif- 
ferential time delay, two approaches are pursued. AR model performance, Doppler es- 
timation and two types of time delay estimation are examined in Chapter V. In the last 
chapter some general conclusions of the work carried out in this thesis are presented, and 
some suggestions for future investigations are given. Computer simulation programs are 


included in Appendices E and F. 


Il. COHERENCE 


A. DEFINITION 
The coefficient of coherence between two wide sense statlonary random processes 
is the normalized cross power spectral density function defined by Wiener [Ref. 1: p. 12] 


as 


Gxy(f) 
a = (27g 
N Gxx/) Cred) 


Where / denotes the fredlucne. (72) 
G,,(f) is the cross power spectrum between x(z) and 3(7) . 
G{f) denotes the aule voy en Speciminy of vi eeciad 
G,,(f) denotes the auto power spectrum of y(z) . 
Despite some confusion in the literature. Wiener intended for the coefficient of co- 


herence to be complex. This is apparent since he discusses both the modulus and the 


argument of the coefficient of coherence. 


B. PROPERTY OF THE COHERENCE? UN GiGi 
The power spectral a matrix Q(f) is positive semi definite. Therefore, for two 


fandonl pracesses « and 3, see that 


G Ga 
det [OW] = det wl) Go) Ea (2.2) 
Gyxff) Gy (/) 


For real processes we have G,,(f) = G},() and thus 


Gf) Gy,(/) = | Gx) | ° 20 ( 


tJ 
Ws 
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where * denotes the conjugate of a complex number. And 
Gxxf/) Gy) Ee | G y | : (2:4) 


Furthermore, G,,(f) and G,(f) are nonnegative, real functions of frequency. When 
G,,.(f) and G,,(/) are strictly positive definite ( 1.¢., G,,(/)G,G) > 0 ). Eq: (2) Canseemmme 
vided by G,,.(G,,(/) without changing the sense of the inequality. 
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This provides as an upper bound 
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and since the magnitude of anv complex number is greater than or equal to zero, we 


have the lower bound 


The magnitude of the coherence function 1s always between Zero and one. It is zero 
if the processes x(t) and y(z) are uncorrelated and it is equal to unity if there exists a 
linear relationship between x(z) and j(z). In order to define the coherence it is necessary 
that the numerator and denominator of Eq. (2.1) are not simultaneously equal to zero. 


Coherence 1s not defined if either auto spectra ts Zero. 


C. COHERENCE ESTIMATION 
Iie 7, denotes the Fast Fourier Transform (FFT) of the / th segment of x(7) at 


frequency f, , then the spectral density estimates are given bv [Ref. 2: p. 22] 





N 
Gxlfi) = a) LMA 2.7) 
i=] 
. 
Cry(fi) a 2) XU Yi (fp) (23) 
a] 
t= 
N 
Cid = x) Vif) | (2.9) 
l=] 
Where o = 
sean. 


eS mUIMoer O1 Seoments, 
T = segment length, and 


/, = sampling frequency. 


Finally the coherence estimation 1s 
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D. COHERENCE OF NARROW BAND SIGNALS WITH DIFFERENTIAL TIME 
DELAY AND DIFFERENTIAL DOPPLER 

The output of the band pass filters (BPF1l and BPF2) in Figure | are denoteqaiag 
Af.) and 1.) respectively. Each term represents the Fourier transform Of the w@omms 
sponding time series evaluated at frequency f, and time 7. For narrow band signals a 
Doppler shift corresponds to a frequency shift. If a signal arrives at the two sensors 
having a differential Doppler shift as well as a differential time delay, then we see that 
a frequency compensation and time delay compensation by the appropriate values tend 
to line up the signals in frequency and time. This is accomplished by using an additional] 
Fourier transform in channel-! of the processor and a trme delay in channel-2 of the 


processor. Nfathematically, this can be expressed as 


v 


» MA Yadide” 
id 


l=] 


rf) oo (2.11) 
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Comparing this with Eq. (2.10), we see that we generalized the coherence concepiws we 
also note that the implementation resembles a correlator, where one of the signals 1s 
frequency compensated and the time delay corresponds to the delay operation of the 


correlator. 
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Figure 1. Coherence estimation block diagram. 


Meier rivtre 1s redrawn as in Figure 2, then it can be mnterpreted as an [FI 


implementation as shown in Figure 3. 
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Figure 2. Coherence estimation block diagram (reinterpretation). 
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Figure 3. Colierence estimation block diagram using the FFT. 


But this FFT has a poor resolution. To obtamn a better resolution, an AR model 


is destrable. This approach ts shown tn Figure 4. 
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Figure 4. Coherence estimation block diagram using an AR model. 


Figure 1 and Figure 4 represent two different schemes to compute the coherence 
function. Differential time delay and Doppler difference can be estimated by using AR 


models in the modeling of the coherence function. 


Wi, AUTOREGRESSIVE (AR) MODELS 


A. AR MODELING 
If the following difference equation 1s satisfied, the resulting structure 1s called an 
AR model of order p. [Ref. 3: p. 177] 


p 
x(1) = — > aex(1 — k) + u(t) (sal) 
k=} 
where x(7) = the signal at instant n, 
u(t) = the white noise driver, and 


a, = the & th coefficient of the AR model. 


A realization of the AR model is illustrated in Figure 5. 
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AR 
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Figure 5. — AR filter of order p. 


B. ADVANTAGES OF THE AR MODEL 

The motivation for parametric modeling of random processes is the ability to obtain 
better spectral estimates based upon the model than estimates produced by classical 
spectral estimation. Both the perrodogram and correlation methods can be used to vield 
Power Spectral Density (PSD) estimates. The unavailable data or autocorrelation se- 
quence (ACS) values outside a given window are assumed to be zero. This kind of an 
unrealistic assumption leads to distortions in the spectral estimate. 


The advantages of the AR approach are 


I. AR spectra tend to have sharp peaks, a desirable feature of high resolution spectral 
estimators. 


2. AR parameter can be obtained as solutions to linear equations. 


C. POWER SPECTRAL DENSITY 


Eq. (G.1) can be rewntten as folliens 


Is 


x(z) =— ) apx(r — k) + ule) 


? “ 
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~ 


% 


= > rade —k) 


(3729 
n=) 
Where 4, = the causal filter impulse response. 
Let us take the Z -transform of Eq(3.2). 
ie — Tie) (3.30 
Feesuittliowisde (os) eas 
p 
> ayx(n — k) = u(r) (3.4) 
k=0 
where a, = 1, and taking the Z -transform gives 
A(z)A(z) = U(z) (3a 


Eq. (3.3) and Eq. (3.5) cai be Solved 10 obtains, 





H(z) = Ale) (3.6) 
and hence 
eee). : 
ie) — Al) (3a7) 


The Z - transform of the output sequence {x(7)} 1s related to the Z - transformation 


of the input random process u{n) by [Ref. 4: p. 56] 
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The AR power spectral density is obtained by substituting z =e?” into Eq. (3.8) 
p F ; } g 


and scaling it by the interval T. 
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P,. = Variance of driving sequence. 


D. BURG’S ALGORITHM 

In practice, the autocorrelation 1s usually not available, so one must make an AR 
spectral estimation based on the available data samples. The simplest procedure to ob- 
tain an AR spectral estimate from data samples would be to produce estimates of the 
autocorrelation sequence from the data. These autocorrelation estimates would be used 
in lieu of the true autocorrelation sequence in the YULE - WALKER equations to vield 
the AR coefficients. However better results are obtained, particulary for short data 
segments, by algorithms that obtain the AR model parameters directly from the data, 
without explicitly estimating the autocorrelation function. 

The Levinson recursive solution to the YULE - WALKER equations relates the 


pth order AR parameters to the p— 1 th order parameters as given by [Ref. 3: p. 211 | 


a(n) = a,_y(n) + kyay_y(p — 7) (3.10) 


For n=] to n=p-—l1. the reflectron coefficients k, can be found bv using the 


known autocorrelation function for lag 0 to p— 1. So we have 


p-i 
a > a, -(A)rg(P = 7) 
n=O 4 
\) - 0?) = (3.11) 


The recursion for the driving white noise variance is given by 
2 
Pp = Ppl — [kp I) (3.12) 


WiNGhe Gaia 10 ie 

But the ACS is not available, hence we can not calculate the reflection coefficients: 
The Burg algorithm provides an estimate of the reflection coefficients which in turm are 
obtained through a least squares criterion. 


The forward linear prediction error is given by 


P 


e(n) = x(n) + > ome — 1) (3.3) 


m=) 


while the backward linear prediction error 1s given by 


ie 


2 (n) = =x(1— p)t+ a, (m).x( n+m-—p) (3.14) 


ae 


x 


= 


Substitution of Eq. (3.10) into Eq. (3.13) and Eq. (3.14) yields the recursive relationships 


et(n) a (1) 7G rae (az — 1) (Sal 
e,(n) = e, (1 —1)+ ke! (n) (3.16) 


At each order p, the arithmetic mean of the forward and backward linear prediction 


error power ( sample prediction error variance) 1s given by 


10 


ec. 2, 1 b 
=| >, leon? +—- dD lest)? (3.17) 
n=p+] n=p+] 


tiswes Pression 1s muninuzed, subject to the recursion given by Eq. (3.15) and Eq. 
(3.16). Thus, p* is a function of single parameter, namely the complex valued reflection 


Seciicient A. Setting the complex derivative of Eq. (3.17) to zero 


ep a ep, 





=a 3.18 
GRe(k,) | ~ lm(k,) >) 
and solving for A, vields 
‘ 
si y ef _y(n)ep_y(x — 1) 
a ee 
H=pt 4 
a (3.19) 
Z 5 
y le_ al’ > » eaeea i) 
n=p+) n=p+] 


The estimation of the reflection coefficient represents the HARMONIC mean be- 
tween the forward and backward partial correlation coefficient. where 


e(72) = e3(7) = x7) and Vs = number of data points. 


E. FINAL PREDICTION ERROR (FPE) CRITERION 

Because the best choice of filter order 1s not generally known a priori, it is usually 
Mecessary 1n practice to postulate several model orders. FPE is a kind of criterion which 
tease provided by Akaike [Ref. 3: p. 230]. This criterion selects the order of the AR 
process so that the average error variance equals the sum of the power in the unpre- 
dictable part of the process and of a quantitv representing the inaccuracies in estimating 


the AR parameters. The FPE for an AR process is defined as 


eee 2) 


wh Pe XT e+) 


where .V is the number of data samples, ys 





pis order of thesititensana 
p, is the estimated white noise variance when using a p th order filter. 


The order p selected is the one for which the FPE ts nunimum. 


IV. DOPPLER AND DIFFERENTIAL TIME DELAY ESTIMATION 


A. DOPPLER ESTIMATION 


Let x(r) and y(t) be the signals received by two sensors 
x(t) = cos{22(f+ a, )z} 
y(t) = cos{(2n(f + a5)2} 


mere) 1s the carer frequency, 
o, .%, are Doppler shifts. and 


FAS the tine waria ble, 


Let f, be the sampling frequency which satisfies the Nyquist theorem. For con- 


venience, in all derivations and simulations a sampling frequency of 64 Hz is used. to- 


gether with band pass filter width of 1 Hz. We assume that 


PAO ou t= 12 yi. 


Which implies that the signals stay in the band pass filter regions of their respective band 


pass filters regardless of anv Doppler shift. Note, these values can be modified to arbi- 


trarv sampling rates and pass band regions. The sampling rate 1s given by 


f, = + 1 > 2Uf+ a) ( r= 1.2 ) 


The phase at instant k and sampling time interval 7 are given by 


2n(f + o,)k 
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Let x(7) denote the sampled analog signal x(r)| _,, then 


x(a + m) = cos{2n(f+ a,) i: + Pm} 


— 
fe 
(42 

Les 


(4.7) 


v(a +m) = cosf2n(f+ a) $ + Pn} (4.8) 
5 


The derivations of Eq. (4.7) and Eq. (4.8) are given in Appendix A. 

These signals. {x(7)} and {y(m)}. are processed to detect the Doppler differeneas 
Let BPFI1 be the band pass frlter centered at f/, and BPF2 be the band pass filter centenem 
at f; (Where f, might equal /, ). VN data points are processed in the band pass filters. to 
generate one output. The inputs to BPFl and BPF2 at time are vectors A, (7)eame 
(2) . respectively. The size of the input vector to the filter is the number of data points 
taken during | second (1.e.. the number of points processed in the filter =—¥ = 7 see 


Input vectors are denoted by 
XA, (1) = (x01). xa + 1)... eG + NV = 1] (4.9) 
Y() = Lv(). 01 + 1)......, yon t+ V— 1] (4.10) 


The filtering is performed using FFTs . where successive FFT outputs are generated at 
the input data rate, The BP? outpul is comjpmeaced: 


To avoid the complexity the following four complex constant variables are defined. 
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A= D cost2alf+ a0) Je” af (4.13) 


as 


oo 


al 


—" 


B= 
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i sin{2n(f+ a) 7 PWN (4.14) 
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At instant #2 the complex output signal of BPF1 can be calculated as follows 
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The complex signal 4,,(/) contains two frequencies with different amplitudes. 
emuncderstand the Character o! -\(/) , it 1s important to evaluate which is the 


dominant term in the above expression. 
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From Appendix C 


As 7s) = ee (4.18) 
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Therefore -——S=-— en 15 the doniunant temo en 
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In sinuhar way, the complex output signal of BPF2 can be calculated as follows 
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The complex signal ¥.(f) contains two frequencies with different amplitudes. 
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Even though }7(f) is similiar in form to_X,(f), 1t 1s not obvious which term is dom- 


inant: 
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From Appendix C 
| A, —jB,| > | A, +/B,| (4.22) 


cal | ' n 
Therefore —~—e'y¢m is the dominant term of Y%(/). 
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Since the two output sequences LY,(/)} and {¥7()} have two sinusoidal components 
each, their product {V,() ¥i()} has four sinusoidal components. The four frequencies 


are 


@, = 22 ———_ ,w), = 2 — 7,03 = , OO, =e 
Is Is 
with @, the sinusoid with the largest amplitude and @, the sinusoid with the smallest 
amplitude. 
The product of the output of the filters is given by 
ore, Ay +/B, 5.9, Ay By eThOn 
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where f(a, 4. .@,.,?,) represents the three low aimiplitude freqnemc; stem 

When using the AR model as described in chapter II], at a high S.\Aee 
SVNR— oo ), a 4th order AR model detects all of the four frequencies given above. 
When the S.VR is low (1.e.. 20 dB ). onlv one donunant sinusoid is detected with fre- 


quenes. 


oO. = > a 
(4.24) 


Since we can detect w, and know the value of f , «, — a, can be estimated by using an 
AR model. 


B. DIFFERENTIAL TIME DELAY ESTIMATION 

To detect and localize a target, it 1s important to find the differential time delay of 
signals arriving at two sensors. Let us assume that signals at the two sensors are as 
given in Figure 6(1.e.. zero differential time delay ) For the remainder of the figures, the 


time axis 1s scaled to be 64 points per second. 
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Figure 6. Receiving signals at sensor 1.2 (SNR= 100, no differential delay). 


We can estimate the differential time delay and differential Doppler using two dif- 
ferent approaches. The first approach uses the cross power spectrum while the second 
approach uses the coherence. 

1. Differential time delay and differential Doppler estimation using the cross power 
spectrum. 

Figure 4 shows the block diagram of the coherence estimator using AR mod- 
eling. In this figure, the output of the FFT can be interpreted as the cross power spec- 
trum. Using this cross power spectrum, we can estimate the differential time delay and 
the differential Doppler. But when we use the highest peak of the cross power spectrum, 
the peak is somtimes not detected at the proper time delay nor at the proper Doppler 
shift. We will show a special case in which the peak of the power spectrum is not de- 


tected at the proper time delay nor at the proper Doppler values. 


a. Special case. 

For our test case the data duration ts six seconds, two seconds of the noisy 
signal and four seconds of noise only. Linear transformation of this data leads to one 
of three types of outputs. The first tvpe represents full information, while the second 
tvpe represents partial information. The third tvpe represents a noise only condition. 
Generally when two signals are hned up in time, the AR model should give the highest 
output power. But this is wrong in some cases. In the full information case, the mag- 
nitude of transformed signal is high and constant. For some reduced information case. 
even though the magnitude of the transformed signal decreases, it still mav have large 
magnitude of spectral components. This phenomenon can be explained as follows. 


Define two functions depending on k 
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Where /p=/f +o. 


A denotes the number of lost data at BPF2 input vector. 


We assume that we know when the BPFI signal starts. If this information 
1s not available. we need to examine the signal at the output of BPF1. to obtain a can- 
didate time frame. Consider the input of AR model as shown in Figure 7. 

The input data size to the AR model is the number of linearly transformed 
data points during a given period (1.e., input size 1s V=/, ). The BPF] output represents 
full information (i.e., each output element is produced from the signal information 


points). Therefore. as shown in the previous section, the dominant output sequence 1s 
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Figure 7. AR model and its driver source. 


(1) O second delay case. We already defined the input vector A,(77) and 
Y(2) to BPF1l and BPF2. This vector can be interpreted as a tine snap shot. Lach 
output of the BPF requires N input points. If we require NV output points from the BPT, 
even using maximum overlap in the processing, we require at least 2\’ — | points at the 
BPF input. Figure 8 shows the 2N input to BPFl and BPF2. Both 2N input points 
contain the signal and some noise. As shown in Figure 8, in the zero delay case, the two 
Input vectors Y,(7) and ¥,(7) (0 < mm < WN ) provide full information. So the lin- 
early transformed output X,(/,) and ¥%(f,) also represent full information. 


The dominant part of the BPF2 output is the sequence 
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Figure 8. Two signals. 0 second delay at the sensors. 


The dominant input sequence to the AR model is given by 
{ky koe et ofan alae ee (4.31) 


where the magnitude 1s £,4, 


(2,  -1 second delay case. Figure 9 shows the 2.V input to BPF1 and 
BPF2. The 2.V input data points of BPF1 contain the signal and noise. During the first 
V points the input to BPF2 contains the signal plus noise while during the second .V 
points onlv noise is present. The input vector X,(”) of BPFi contains full information, 
but input vector Y¥.(7) of BPF2 does not contain full information. In the },(7) case, 
every element of ¥,(7) is full information. But when A is not equal to zero then Y,(7) 


consists of \ — & information elements and & noise elements. 
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Figure 9. Two signals, -1 second delay at the sensors. 


The output from BPF2 represents partial information while the out- 


put from BPF] represents full information. The BPF2 output sequence is given by 


nay, Tp By rty =F) 4B 


2 — Oe 4 an = Oh re) eee, = 1 i} (4.32) 


The derivation of Eq. (4.32) is given in Appendix B. 
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Assuming that “+ e~/y*k 16 the dominant term, which is a valid assumption for our 
fest case with f =6482z, f=23Hz, «,=0.4Hz and o,=0.001 Hz, then the input se- 
quence to the AR model 1s donunated by 
Aj. B, , 
Oty ed i - : 
{ky a oP) = 01.2.0N 1 (4.33) 
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with k, is defined earher ( Eq. (4.28) ). 


The magnitude of dominant input sequence is obtained by examining 
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When we compare the magnitude of the sequence in Eq. (4.27) and the magnitude of the 





dominant term in Eq. (4.32), the magnitude of |4,| > |,4, —j,8,]/2 where 
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We note that Eq. (4.35) has two frequencies (1.e., 0 and > 2m ) 
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The magnitudes of both terms are :. 
| —evx| 


: ] : 
Now, note that as long as |x] < G° then the magnitudes of both terms in Eq. (4.35) 
) 


are larger than |4,| 


Due to the svmmetrv a positive delav results in a simular reponse to a negative delay, and 
hence Eq. (4.35) holds also for delay of +1 second. Equation (4.35) can be interpreted 


as having large responses at shifted Doppler values at a delay of | second and a delay 


of —1 second. Proper delay can not be established nor can the Doppler value be estab- 
lished. This cross power spectrum algorithm will not give good information about either 
Doppler difference or time delay. 

b. Nlodified cross power spectrin. 

To detect the signal, input signals to the AR model should be preprocessed 
as shown Figure 10. If the BPF2 output signal 1s magnitude normalized, then this nor- 
malized signal retains good phase information. When we normalize the BPF2 output, 
one of two conditions can occur. The first condition ( 1.e., synchronized ) leads to the 
magnitude normalization of Eq. (4.29) and provides accurate frequency and delay esti- 
mation. The second condition occurrs when the signal in channel-2 is not lined up in 
time with the signal in channel-] ( 1e., during the delay search ), Under this condition 
smaller peaks in the tme Doppler plane appear at incorrect values of time delay and 
Doppler. Above 70 dB SVR, the correct peak 1s the dominant one and allows proper 
estimation. Information from contour plots can be used at values of SVR between 70 
ama 10 dB ( see Appendix D ). 
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Figure 10. Modified cross power spectrum block diagram. 


2. Differential time delay and differential Doppler estimation using the coherence. 
Another way detecting the signal is normalization of the AR output with the 
squared root of the product of P,,(/) ( auto spectrum of x ) and P,,(f) (auto spectrum of 


y ). This can be done using two additional AR models as shown in Figure 11. 
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Figure 11. AR model of coherence. 


Using the AR approach, p,, , Pry + Py , aNd AR coefficients can be Ob: 
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where P_,(f) is the power spectrum of A, , 
Pf) is the power spectrum of }7_,, 
emits tne power spectrum of 1,27; , 
A,,(f) is the AR power spectral density of .4, , 
A,,(f) is the AR power spectral density of Y, . 
A,,(f) is the AR power spectral density of 4,%)., cross term. 
Peete Ine driving Moise vamance Of A, : 
Py, 18 the driving noise variance of 17, . 
Psy IS the driving noise variance of A,¥7_,. and 
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We obtain a coherence estimate from the three power spectra, by assuming that 
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This approach provides good information about the differential time delay and 
eierential Doppler down to S.\Rs of 20 GB. 
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This chapter presents graphical results obtained by applving the analytical results 
of the previous chapters to specific examples and carrying out the required computations 


on the computer. 


A. AR MODEL 

This section provides the AR model performance tests. As we discussed in chapter 
Ili, Burg’s algorithm is apphed to this AR model. The FPE criterion is use to find the 
AR model order. Figure 12 1s the power spectrum of two test signals sin(w,/) and 

os(w,/m). Figure 13 shows the power spectrum of two other test signals sin(@,77) and 
os(w,/”). Both sets of test signals in these figures have an S\R of 20 dB. The free 
quenev f, 18 13.45 Hz and fi is 23.45 Hz and Tis 7 second. 

Theoretically, the power spectrum of these signals have two impulse functions at fi 
and —/f, ( f=1,2 ). Figures 12 and 13 show experimental results agreei@ witli 
theoretical results. Figures 12 and 13 indicate which frequency has high power Duties 
provide no phase information. 

pee: AR model is sensitive to frequency but 1s not sensitive to the phase. For anv_ 
selection of frequency f which satisfies the Ny Guist theorem, the AR model provides 


good frequency estimation provided the S Wk as * suflicient large. 
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Figure 12. 
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AR model performance test | (SNR = 20 dB). 
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Figure 13. AR model performance test 2 (SNR= 20 dB). 


30 


0 20 


FREQUENCY(HZ) 


=20 


B. DOPPLER ESTIMATION 
In section A of chapter IV, we discussed five items. 


1. BPF1 and BPF2 outputs have two sinusoidal components each. 


2. For BPF1. the dominant frequency is designated as w, and for BPF2, the dominant 
frequency is designated as w, where 
— 2+ 1) 
O, = _— 
8 —2u/ + %) 
Oy = ——- 
3. The power spectrum of cross term of BPF1l and BPF2 has four components which 


can be detected in very high S.VR (.e. SNR = 100 dB). 


pemeveien the S.\Ak is low (i.e. S\R = 20 GB). only one dominant frequency compo- 
nent is detected. 


5. The dominant frequency of the power spectrum P 
ential Doppler frequency. 


waxy) Corresponds to the differ- 

ficmre 14 shows the power spectrum of the BPF1 output when fis 13 H:. o, 1s 
0.45 Hz, and fis 64 Hz. Figure 15 shows the power spectrum of the BPF2 output 
eememe, 15 13 Fz. «, 1s —O0.45 Hz. and f, is 64 Hz. As expected both figures show two 
spectral lines when the S.VR is 100 dB. The dominant frequencies are located at f equal 
free iz tor BPF1, and at f equal —12.55 Hz for BPF2. Figure 16 shows the power 
spectrum of the cross term of BPF1 and BPF2. This figure shows four spectral lines at 
an S.VR of 100 dB. while only one dominant frequency is detected at an S.VR of 20 dB. 
Domunant frequencies are alwavs located between —1 Hz and 1 Hz. Figure 17 shows 
a subplot of Figure 16 for frequency between —1 Hz and 1 Hz. This figure shows the 
dominant frequency f=.9 Hz and the smaller spectral component at f= —.9 Hz for an 
S.VR of 100 dB. Only one donunant frequency can be detected at fequal 0.9 Hz for an 
SNR of 20 dB. Doppler difference from the dominant spectral component corresponds 
mortne true Doppler difference «, — a, which ts 0.9 Hz 

Figure 18 shows a comparison for SVR of 20 dB and 10 dB. Figure 19 is subplot 
of Figure 18 for—1 Hz</f<1 Hz. From this figure, when SNR = 10 dB, the dominant 
frequency is shifted a little bit from the true Doppler difference. We have shown the five 
issues as We have discussed in the previous chapter. In summary, for S\Rs_ greater than 
or equal to 20 dB, the power spectrum of the cross term provides good Doppler esti- 
mation. 
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Figure {4. 
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Power spectrum of the BPF Ii output. 
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Figure 15. Power spectrum of the BPF2 output. 
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Figure 16. Cross power 
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spectrum of BPFi and BPF2. 
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Figure 17. Subplot of the cross power spectrum. 
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Figure 18. . Comparison the Doppler estimation at two different SNR’s. 
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Figure 19. 
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Subplot of the Doppler estimation. 
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C. DIFFERENTIAL TIME DELAY ESTIMATION 

In Chapter IV section B.l.a we discussed the problems associated with using un- 
normalized quantities and suggested two possible approaches to solve the problem. In 
this section we will Shiow: that for tine caserO1 aon = A Hz , the power spectruniomee 
nondelaved signal is smaller than the power spectrum of a signal delayed by one second 
( i.e.. see discussion following Eq. (4.35) ). 

To estimate the time delay. we suggest the two different approaches. The first ap- 
proach is the normalization of the BPF2 output. The second approach normalizes the 
cross term P,,,,(/f) with the power spectrum of BPF1I P,,(/) and the power spectrum of 
BPF2 P,,() through the equation defined by Eq. (4.40). 


The following figures are simulations for f= 23 Hz, o, =0.23 Hz, #,=—0.02 Hz. 


l 
i 4 


P.,(f) while Figure 21 shows the contour plot of P 





second, delay of 0 and an S.VR of 100 dB. Figure 20 shows the surface plot of 
wy) Figure 22 shows the eyame 
section plot of Figure 20 with respect to the DELAY axis. This figure shows that the 
power spectrum has a spurious peak at delay of one second. At the proper delay (1.e.. 
Q second delay) P,,,,(/) has a relatively smatl value. 

Pil) is affected by two terms. namely p, and 


. Lae 
maXimium power spectrum term of the transfer function and Figure 24 shows the driving 


. Figure 23 showsaune 


noise variance. At a delay of + 1 second and of 0 second the transfer function shows a 


peak while the variance of driving noise has small value. In this case, P,,,,(/). the product 


X¥XY 


Of thesiomnce. at 0 second delay has a smaller value than that at anv other 


asl 
PaO: | 
delav location. Using the scheme shown in Figure 10, we normahze the BPF2 output 
and present the result in Figure 25. We can detect the proper time delay and the prea 
Doppler difference. Figures 26 and 27 show the corresponding contour plot and the 
power spectrum of the normalized cross term. For anv value of time delay and Doppler 
difference, this approach provides good information about differential Doppler and dif- 
ferential time delay provided the S.VA is greater than or equal to /0 dB. 

Results of the second approach ( Figure 11 ) are demonstrated in Figures 28 and 
29. This approach also gives good information about differential time delay and differ- 
ential Doppler. This approach provides good time delay and Doppler estimation for 
5.) hecieater than or equal te 20 al. 

If we use the contour shape, we can also estimate the differential time delay and 
differential Doppler for S.VR of less than 20 dB. Differential time delay and differential 
Doppler estimation can be extended down to S$.VRs of about 0 dB when using the con- 


tour of the coherence surface. Figure 30 is a contour plot of using this second approach 
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fee econerence appioacn ) for &,—90.23 Hz.a0,——90.02 Hz, Delay=0, and an S.\VR of 
O dB. The contour shape of the second approach has the shape of the letter A. When 
the S.VR is greater or equal to 0 dB, the proper time delay and Doppler difference seems 
to be located at the cross over point of the A. When the S.VR ts less than 0 dB, the 
contour shape does not shows an X clearlv. In the coherence approach. Figure 3] 
provides the contour plot, with channel-1 containing onlv noise and channel-2 contain- 
ing signal plus noise ( $S.VR=0dB ). Figure 32 shows the contour plot, when channel-2 
contains only noise and channel-1 contains signal plus noise ( SVR=0O dB ). 
Figure 33 shows the contour plot, when both channels contain noise only. The above 
three figures do not show an X clearly, hence do not allow signal detection for esti- 
mation of anv parameter. Using this information. we can distinguish the signal combi- 
nations from the noise combinations. 

When we estimate the time delav and Doppler difference using the contour plot. the 
modified cross power spectral approach has some problems as discussed in 


Appendix D. 
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Surface plot of the power spectrum. 


Figure 20. 
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Figure 21. 
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Contour plot of the power spectrum. 
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Figure 22. 


Nlaximum power spectrum. 
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Figure 23. 
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Maximum power spectrum of the transfer function. 
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Variance of the driving noise. 
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Figure 25. Maximum modified cross power spectrum of the transfer function. 
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Power spectrum of the transfer function (surface). 


Figure 26. 
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Figure 27. 





OS 0 OS— Oot 
SIXV AV145G 


Power spectrum of the transfer function (contour). 
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Surface plot of the coherence. 


Figure 28. 
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Figure 29. 
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Contour plot of the coherence. 
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Figure 30. 
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Estimation using the contour shape (low SNR). 
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Figure 31. Contour plot (noise in channel-1 only). 
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Contour plot (nbise in channel-2 only). 
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Figure 33. 
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Contour plot (noise only). 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


In this thesis, the differential time delay and differential Doppler are estimated using 
two different approaches. The first approach uses a modified cross power spectrum al- 
gorithm while the second approach uses a coherence algorithm. In both cases 
autoregressive modeling of the required cross spectral component is used. The differen- 
tial time delav and Doppler difference can be estimated using a threshold at high S\ Rs 
and a contour plot allow some 

At high S.VAs ( te., SVR = 70 dB ), the first approach locates the dominamiumems 
at the proper trme delay and the proper Doppler. The second approach utilizes two 
additional AR models and two additional FFTs to obtain the autopower spectrum of 
channel-! and channel-2. At moderate SNRs (1.e., SNR = 20 GB ), the second apprenes 
has the highest coherence peak at the proper time delay and Doppler frequency. When 
the S.\A = 0 dB, the contour plot of the coherence has the shape of the letter Xouaie 
differential ume delay and Doppler difference can be located at the cross over point of 
the A. In high S.VAs we can use the highest peak to detect the time delay and Doppies 
frequency using either approach. When we use the contour plot, the information about 
tume delay and Doppler depends on the location and the form of the cross over point 
Ol tie. 

We can estimate the differential time delay and differential Doppler using the AR 
modeling of coherence. But we can not get the numerically correct value of the coher- 


ence coefltcient. Ihe following three pormts are recommended tOmiutiinresrudn. 


5 
° 





[. To get a cross power spectrum, we assumed inate). (is eee This was 
required since We can not get a cross power spectrum using an AR model. Im- 
provement of the cross power spectrum estimates using an AR model can lead to 
improvement of the algorithm. 


2. Our first approach is not normalized by the auto power spectra. If we find the 
corresponding auto power spectra, We can also obtain a coherence coefficient esti- 
mate. 


To get a better information at a low SNR, we need to study the charactenstiesma. 
the contour plots. 
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APPENDIX A. PHASE DERIVATION 


Define the phase at instant &. 


2m + Ok 


xP = ie 


2n{f + 04)K 


y? k= i 
5 


Let x(z) denote the sampled analog signal x(r) | ,.,7 


then 
x(7) = cos{2a(f+ x) 7; 
= cos{2n(f+ ¢,)—=+ bo} 
x(a + 1) = cos{2z(f+ o,) = a ! } 

a 
= cos{2z(f+ a,) i a 
= cos{2n(f+ a) $ ss fey 

<i 2) = COS{ 2alf + 0) a a 5 } 

Ualeae 2) ie 
= cos{2: aa one 
= cos{2n(f+ a) $ + xP} 

Peer 
x(n + m) = cos{2z(f+ a) aa } 
= cos{2n(f+ a,)— 7 i 
= cos{2n(f+ o,) = oe, 
Simuhianiv 
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(4.7) 


y(n) = cos{2a(f+ 2) 


y(n + 1) = cos{2a(f+ a) 


= cos{2x(f+ a.) —+ 
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= cos{2a({f+ a>) 7 + y@o} 


n+] 
a 
Js 





7 


= cos{2z(f+ a>) F + >} 


y(t + nt) = cosflarff+ a, 


Mo Pl 


eae 


22lf+ O15) rh 
1 er 


Is 


Ss 


Dal f + 01>) 


(a.6) 


(4.8) 


APPENDIX B. OUTPUT OF BPF2 


Define two functions depending on k 
N=1—k 


ety = > cos(2nfy : a (4.25 
n= 


(7 
ee” 


N-\-k 


B= >. sinQnf : Hy Ph (4.26) 


n=0 . 


where f.=/f+o,. 


FULL INFORMATION CASE (N POINTS ) 


. points are data points at BPF2 input. 


N=! 


KW = ) cosa T+ bye 


ow 
i 


a Def 
{ cost 2h COS Oy SIN(z a F) sin Ope x 
n=0 5 

N=] a 


ofr er: , ny, pet 
= = 05 Bm, > cost (2nfy— Je V— sin ym > sin(2nf; f. Je (6.1) 
n=0 


f= 


— oA, COs = 0B, sin - 


hom 4 chin gin gyn 
=i. — 0B, oF 


= oy = B, alyom i 


hn _ 


oy SD y ois o 
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B. PARTIAL INFORMATION CASE ( N-1 POINTS ) 
.V— 1 points are data points at BPF2 input. | 


N=-1-1 
- H apt 
a me » cos(2z 2 ie ae yPme y ‘ 
n=0 5 
N=-1-1 
as af 
= { cos(2zf5 t ) cos .0,, — sin(2xf, 7 —) sin Oe oa 
n=0 
N=-1l=-] ‘ a . 
} Qa : . ri —/) aS 
= COS. Or y ee N ashi) » sin(2z Fe pam (6.2) 
n=0 2 r=) J 


= (71, COS, — 7, or 


_ given eu o fyOm Be mm ot yOm 
= ey > as dj 
es ey 
Byes : 
1 dy, TH; B. ia ] i, Jj B, jo 
= — A ee 


fm — 


C. PARTIAL INFORNIATION CASE (N-K POINTS ) 


\—& points are data points at BPF2 input. 


N-l—z& 


” y ae H \ jiaf 
Pal) = ee cos(2zf; va = Pome ‘ 


n= J 


, so FL ; yl : I-af— 
are OS(2 75 re €0S,.0,, — Sinkeag a ) sin yPmbe 
_ 5 


n= 


N-1—k N=-1—k 
Se _p ny pet » j= 9 
= COS Py, is cos{2zf; . a sin oe Sil(2 7) oan FE V (6.3) 
n= me 20 
> gly COS Om) By sin 
ghsOm =e 5S yO m cl yOm — eo vOm 
= 4- =|, ae 


oo 
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APPENDIX C. AMPLITUDE COMPARISON 


Ecc yand «, be the complex amplitude: 


Gi lee (orl) 


1 
i oe a 


| wo ee 
(95 — SL (c.2) 


where |a&| <0.5, and N > 2f+ 1 > 2(f+ a). The numerators of |C,| and |C,| are same, 
because the magnitude of the complex conjugate is the same. If we show that 
lt — =| < | 1 — ete) | , then we can say |C,| > |C,|. As shown in Figure 34 
Ji — eH] = |All] and [1 — ee y| = [BI . 


0; =—eiroc 


Oo= -271(2f+a)/N 





Figure 34.  Mfagnitudes of two complex number. 


2n(2f+ a) . 2n(N — 2f— a) 7 2n(1 — a) 


The smallest angle Z/OB is O : 
> nN N N 


For positive f 


2n(2f+ a) 


LIOB = 


ao (cr5)) 
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Since | — x] > [a]. 





27 | — oe [oes 
A es 








From Eq. (c.3) and Eq. (c.4) |A/| < | B/| , therefore 


Vere eal 
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APPENDIX D. CONTOUR PLOTS OF MODIFIED CROSS POWER 
SPECTRUM 


mojo cs) ereater than or equal to 10 db the contour plot of the modified cross 
power spectrum algorithm has the shape of a boomerang. The shape of the boomerang 
depends on g,. Figure 35 and Figure 36 show the contour plot of the modified cross 
power spectrum. The two figures show the boomerang shape and their different di- 
feemons,. FOr positive o, { case | ) and negative o, ( case 2 ) the opening of the 
boomerang 1s toward the left and right. respectively. In the both cases, we can estimate 
the time delay and Doppler frequency using these contour plots. The proper time delav 
and Doppler difference can be detected by locating the point of symmetry of the 
Boomerane. Whe Figure 3/ shows a contour plot when o, 1s —0.02 Hz. Since a, 1s verv 
smal], the contour plot tooks like a straight line. In this case, the differential time delav 
and Doppler difference can still be estimated bv locating the point of symmetry of the 
Peemerane { 1.¢., center point }. Figure 38 shows the contour plot. when channel-1 
contains noise only and channel-2 contains signal plus noise (1e.. SVR = 10dB ). Fig- 
ure 39 shows the contour plot. when channel-2 contains noise only and channel-1 con- 
tains signal and noise ( 1.e.. SNR=10 dB ). Figure 40 shows the contour plot. when 
both channels contain noise only. Using the modified cross power spectrum technique. 
Miewthree types of noise only set ups appear as similiar plots ( see Figures 38. 39 and 
40 ). Figure 38 seems to indicate the presence of signals in both channels. while Figures 
Byeand 30 appear more like noise. Hence. contour plots of the modified cross power 
spectrum do not allow clear identification. In this case. it will be difficult to distinguish 
between the noise only situation and a small @, situation. 

item Use (Ne medined cross power spectrum algorithm at low S.\Rs ( i.e.. 


peu /U ), we need tO investigate the contour plot some more. 
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Contour plot of the modified cross power spectrum (case 1). 
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Figure 36. 
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Contour plot of the modified cross power spectrum (case 2). 
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Contour plot of the modified cross power spectrum (case 3). 
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Figure 38. 


SNR=10 CHANNEL1=NOISE 


DOPPLER AXIS 





O01 OS 0 OS— co. CC 
SIXV AV13C 


Contour plot of the modified cross power spectrum (noise in channel- I 


only). 
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Fiotines3?. Contonr plot of the modified cross power spectrum (noise in channel-2 


anly). 
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Contour plot of the modified cross power spectrum (noise only). 
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APPENDIX E. MODIFIED CROSS POWER SPECTRUM PROGRAM 


waeuleniontsnicuen'sn'esltsntsalsnionsn'au'se'se'sa'snivelen' sn’ ow a wlam’, owlontan' onsen entantou'en’ suloute So oso oa tele s‘ ese Se ver'esiesle we s'este a sje ex's s'e ese ele ve se serie verte slew slevles! aa 
OR FL EL EN 4D ED EX ER FRED ER ER CR ER FR EN EN FR ER FR EV ERX FLT ERN ER ER EN EU ED EL ED EXER EX EX FY 3 fy ey eh wee 


ec 

c ss 
C main program of approach one. we 
e this program share subroutine with a 
c the main program approach two ¥ 
C w 
CTE TET ETE TE ETT EEE TE TENS ETE ETE TE TET TE TEE TE TOTS TS TET TE TOTS TOTS TCT TE TST TESTES TET TEETER TER EEE ERE EERE TE 


COMPLEX X(- -128: 255) ,Y(-128: 255), ade 2500), B(0: 2500) 
COMPLEX BPF1(-128: 191) ,BPF2(-128: 191) ,BPF(0: 100) 
INTEGER SNR,BINDEL,DELAY 

DATA X,Y,A,B/5770%(0. ,0. )/,BPF1, BPF2, BPF/741*(0. ,0. )/ 
OPEN(UNIT=7,FILE='FILENAME FILETYPE A‘) 

PI=ACOS(-1. ) 


Ct ttt ete eee me meee em mw mew mmm mee eww wee ee ewe ee ee ee eee eee eee ee ee ee ee ee eee ve 
e > 
c input part We 
C 7 
e fs - sampling frequency e 
E f - carrier frequency 
e alphal - doppler effect at sensorl ve 
e alpha2 - doppler effect at sensor2 * 
C delay - signal receiving time difference we 
ce between two sensors * 
E SNR - signal to noise ratio ¥ 
c iseed - noise generator seed(odd number) ¥ 
C n - number of data during one second * 
Cc se 
Ctr tt tert rrr eter r tre terete ese errs ester eres teres reer er eee eeeeeereeeeeeece=-= 8 

FS=64 

F=23. 

ALPHA1=. 23 

AEP EA 2 =,.02 

DELAY=0 

SNR=100 

N=04 

Ta2 ler Es 

ISEED=13 
Co ccc eee een nee nee ee ne ne eee ee ee ee ee ee ee eee eee ee eee ee ee ee eee ee ee eee eee 3 
Cc wt 
Cc input data sampling at two sensors * 
C “ 
Coc ccc ce ee ee eee ee nw ee ee ee ee ee en ee eee eee ee eee eee eee we 


CALL SIGNAL(X,F+ALPHA1,FS,0,SNR,ISEED,0,0) 

CALL SIGNAL(Y,F+ALPHA2 ,FS,DELAY,SNR,ISEED,0,0) 
DO 100 K=-128,191 

DO 50 J=K,K+63 
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Cc =a 
e linear transformation using band pass filter ws 
Cc ay 
Cotte en ne ee nw nn a ee ee ee ee eee eee eee eee “ 
OMEGA=T*FLOAT( J-K) 
BPF1(K)=BPF1CK)+X(J)*CMPLX( COS( OMEGA) , -SIN(C OMEGA) ) 
50 BPF2(K)=BPF2(K)+Y(CJ)*CMPLX( COS( OMEGA) , SINCOMEGA) ) 
Crete ttee woes eee eee eee ee eee eee ee ee ese ew ee em et ee etm ee eee te ester eee eee rere-e- ae 
Cc ve 
c normalization of BPF2 % 
Cc ve 
Cot ct cee meme me em we ne ee ee ee ne meee eee ee ee ee ee ee eee ee eee ee ee eee ee eee ee ee v 
100 BPF2(K)=BPF2(K) /CABSCBPF2(K) ) 
DOr S00 BIND El—-128 5127 
DO 200 K=0,N 
Cote eee mem mee meme ee eee ee eee eee ene nen ee ee ee ee eee eee ee ee ee ee ee ee ee vv 
ec 
€ input of AR model 
Cc 2° 
ec eerrenenre2ce2e2 ene envZ w@ewewe2ekwe2e2e ee e2 2 e2e2r2ze2ee2 ez 2w@eeerewewe2een2 ee e2 e2 2e2 2@eee2enee2 ee e2eeee2e2 @ @ @ @ @ @= i; 
200 BPF(CK)=BPF1(K)*BPF2(K+BINDEL) 


CALL CLEAR(A,B, 2500) 
CALL BURGAR( BPF ,N+1,A,IP,VAR) 
CALL CHANGEFFT(B,A,11,MAX) 
DO 250 I=-32,-1 
250 WRITE(7,997)BINDEL,FLOAT(1)/32. ,CABS(1. /B(2048+I1 ))*VART 
* ,CABS(1. /B( 2048+1)) 
DO 300 I=0,32 
300 WRITE(7,997)BINDEL,FLOAT(1I)/32. ,CABS(1. /B( 1) )**VART 
* CABS(1. /B(1)) 
CLOSE(7) 
STOP 


oo] ee eid ete 2 bo. 6,26,613.6,2X,E13. 6) 
END 
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APPENDIX F. COHERENCE PROGRANI 


vs 

f h tw ¥ 
main program of approach two sf 
7 


fe cleclecteslocteslonte cleclo stestes'e sleves's s'e sleste s'e s‘e ste s'e ste s'e sleste steve ste estes s‘s sieves's s': s‘e sies'eve s'e Jostesteste sles'esle s'e stes'e es'e 


COMPLEX X(- =1265 255) We 26r Zoe ACO: 3500), BC 0: 2500) 


COMPLEX BPF1(-128: 191) ,BPF2(-128: 191) ,BPF(0: 100) 
COMPLEX AUTOX(0: 100) , AUTOY(0: 100), SUMXY 

INTEGER SNR,BINDEL, DELAY 

DATA X,Y,4,B/5770*(0. ,0. )/, BPF1,BPF2,BPF/741*(0. ,0. )/ 
DATA AUTOX, AUTOY/202"*(0. ,0. )/ 

OPEN( UNIT=7 ,FILE='FILENAME FILETYPE A’) 

PI=ACOS(-1. ) 


input part 


fs - sampling frequency 

fs - carrier frequency 

alphal - doppler effect at sensorl 

alpha2 - doppler effect at sensor2 

delay - signal receiving time difference 
between two sensors 

SNR - signal to noise ratio 

iseed - noise generator seed(odd number) 

n - number of data during one second 


ola 
ry 


@ee@eeeee2@8 @ @ @ ew ewee298 2e28e2 829 #2 ew8 BF BF FB Fwewew FBV’"FVZBEZTHXZTZ7, SFeg EF FT ZF" FF "BPeewreewe wesw wZwZ’_sFsegxiOe #»»e#*@®Weweweeweewewaeswdseoerwteeweswaeswesw sw Bs se @ @ «¢t 


ALPHA1=. 23 
ALPHA2=-. 02 
SONR=20 

F=23: 

N=64 

FS=64. 

Tae ies 
ISEED=23 
DELAY=0 


CALL SIGNAL(X,F+ALPHA1,FS,0,SNR,ISEED,0,0) 

CALL SIGNAL(Y,F+ALPHA2 ,FS,DELAY,SNR,ISEED,1,0) 
DO 100 K=-128,191 

DO 50 J=K,K+63 


@eeeeee2 @ es ess @Beweeeeeseeor"2 @2eeweeee@eoeeeeewee@eeweeeees @eeese# ewe@eeeoeeee2@2 @2¢@2° 2 22 28 82 8 @ ¢ 


efeeeweeeeeeeeeoeee@ #28 Bs @B sw BB" ws BF ew ww wT FTF BT BW sZTZ’_HX Fee »*»e*Weweooeooe#se##w#eq@eoeweertaewewesw FT FBT FT FT FT "X V" FB " ss |" FF @ |] @s |] ss 
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250 


300 


297 


- 
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OMEGA=T*F LOAT( J-K) 
BPF1(K)=BPF1(K)+X(J)*CMPLX(COS( OMEGA) , -SIN(C OMEGA) ) 
BPF2(K)=BPF2(K)+Y(J)*CMPLX(COS( OMEGA) ,SIN( OMEGA) ) 


CONTINUE 
SUMX=0. 
DO 150 K=0,N 
ot ae ee nan an an nw on one nanan an moe ena waa eee eee sk 
find AR order at sensorl oe 
x 


AUTOX( K)=BPF1(K) 
CALL CLEAR(A,B, 2500) 
CALL BURGAR( AUTOX,N+1,A, IP, VARX) 
CALL CHANGEFFT(B,A,11,MAX) 
XMAX=CABS( 1. /B(MAX)) 
DO 300 BINDEL=-128,127 
DO 200 K=0,N 


eee eweweweeweeweweweeewee2we2ne2nnrne2n ewww e2nwnn2newewne2enre2ee2eerwe2eeweewnee2we2eewwrwrewwrwewe2wenwne2ee2eewreeweeeweweeew ew @ ww ow = 3s 


AUTOY(K)=BPF2( K+BINDEL) 
Ber K)—BPPICK)*(BPP2(KEBENDEL)) 


ae 
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oe 
@e@gmqge@ge@®e@a@dq@e@eq@eq@q@e@q@e@q@eq@e@eq@e@®e@e@eq@eq@eqe@eq@eq@eq@mde@eq@e@eqe@eqe@eq@eq@eqe@e@aq@dee@eq@q@q@e@eq@eaq@aq@de@eq@eq@e@eqeeq@ge@eqe@eq@eqe@e@eqgde@eqg@@ si, 


CALL CLEAR(A,B,2500) 

CALL BURG( AUTOY ,N+1,A,IP, VARY) 

CALL CHANGEFFT(B,A,11,MAX) 
YMAX=1. /CABS( B(MAX)) 


ae 
afaoeeweeeee =# ese 2 ese 2B ee ses ese eee ewewewewese® @eeeeaeeweenwrwreg@meeew@ewrerwewwewgewewrereewrteeeeeaeeeeweeeE2 & ws = = = ss. 


CALL CLEAR(A,B, 2500) 
CALL BURG(BPF,N+1,A, IP, VARXY) 
CALL CHANGEFFT(B,A,11,MAX) 
XSMAX=( XMAX**YMAX ) "2°" VARX*"*VARY/VARXY" fs 
DO 250 I=-32,-1 
WRITE(7 ,997)BINDEL, FLOAT(1)/32. ,CABS(1. /B(2048+1) ) /SQRT( XSMAX)*8. 
* CABS(1. /B( 204841) ) 
DO 300 I=0,32 
WRITE(7,997)BINDEL,FLOAT(1)/32. ,CABS(1. /B(1)) /SQRT( XSMAX)*°8. 
* ,CABS(1. /B(I)) 
CLOSE(7) 
STOP 


FORMAT(1X,14,2X,E13. 6,2X,E13. 6, 2X,E13. 6) 
END 
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C2 OO) CG) 3 - GS ae 


100 


C2) CC Gu) a). CG) 


SIGNAL GENERATER 


INPUL 
N - NUMBER OF DATA POINTS 
F1,F2 - FREQUENCY ©F 157,2ND sie) 
AMP1,AMP2-FREQUENCY OF 18T,2ND SIGNAL 
OUDECE 


ACN) = SLGANL 


SUBROUTINE SIGNAL(B,F,FS,D,SNR,ISEED, ID,NOISE) 
COMPLEX B(-128: 255) 
REAL A(0: 1) 
INTEGER SNR,D 
PI=ACOS(-1. ) 
DATA A/O. ,0. / 
SIGMA=1. 
A(0)=SQRT( 2. **10. 2*(SNR/10)) 
IF(NOISE. EQ. 1) A(0)=0. 
DO 100 I=-128,255 
M=(I-D)/128 
IF(I. LT. D. OR. M. NE. 0) M=1 
T=AMOD( F**FLOAT( I-D) ,FS) 
CALL GAUSS(ISEED,SIGMA,0. , RANDOM) 
IF(ID. EQ. 0) THEN 
X=COS(2. *PI*T/FS)*A(M)+RANDOM 
ELSE 
K=SIN( 2. *PI*T/FS)*A(M)+RANDOM 
ENDIF 
B( I)=CMPLX(X,0. ) 
RETURN 
END 


BURG ALGORITHM 


PSP ad 
N - NUMBER OF DATA POINTS 
X - INPUT SIGNEE 
OUTPUT 
IP - ORDER -OPSAR 
ACO: IP)="AR COEPPUGIENES 
VAR - DRIVING NOISE VARIATION 


SUBROUTINE BURGAR(X,N,A,IP,VAR) 

COMPLEX X(N),A(0:N),EFK(500) ,EBK(500) 

COMPLEX EFK1(500),EBK1(500) ,AA( 20,20) ,SUMN, SUMD 
REAL RHO(0: 80) ,FPE(0: 80) 

INTEGER START 

RHO(0)=0. 

FPE(0)=1. E64 

ACO)=CMENA( 1-60.) 

IP=1 
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100 


Crea G)-C) Cy Cle)" Cs Caer ea. co 


START=1 
DO 10 I=1,N 
RHO( 0 )=RHO(0)+CABS(X( 1) )***2/FLOAT(N) 
on Om i=2eN 
EFK1(1)=X(I) 
EBK1(I-1)=X(I-1) 
LOOP 
K=IP 
SUMN=CMPLX(0. ,0. ) 
SUMD=CMPLX(0. ,0. ) 
DO 30 I=K+1,N 
SUMN=SUMN+EFK1( 1)**CONJG( EBK1(I-1)) 
SUMD=SUMD+CABS( EFK1( I )***2)+CABS(EBK1( I-1)**2) 
SUMD=SUMD+CABS( EFK1(I) )**2+CABS(EBK1(I-1) )?*2 
AA( K,K)=-2. *SUMN/SUMD 
TEMP=1. -CABS( AA(K, K) )***2 
IF(TEMP. LE. 0.) TEMP=1. E-10 
RHO( K)=TEMP**RHO(K-1) 
IF(K. GT. 1) THEN 
DO 40 J=1,K-1 
AA( J,K)=AA(J,K-1)+AA(K,K)*CONJG( AA(K-J,K-1)) 
ENDIF 
DO 60 I=K+2,N 
EFK(1)=EFK1( 1)+AA(K,K)*EBK1(I-1) 
EBK( I-1)=EBK1( I-2)+CONJG( AA( K,K) )*EFK1(1I-1) 
DO 70 I=K+2,N 
EFK1(1)=EFK(I) 
EBK1( I-1)=EBK(I-1) 
IF(N-K. EQ. 1) THEN 
FPE(K)=FPE(K-1)+1. 
ELSE 
FPE(K)=RHO( K)*FLOAT(N+1+K) /FLOAT(N-1-K) 
ENDIF 
IP=IP+1 
UNTIL( FPE(K). GT. FPE(K-1). AND. K. GT. START) 
IP=K-1 
DO 100 I=1,IP 
A(I)=AA(I, IP) 
VAR=RHO( IP) 
RETURN 
END 


BURG ALGORITHM 


INPUT 
N - NUMBER OF DATA POINTS 
X - INPUT SIGNAL 
re - ORDER OF AR 

OUTPUT 


PCO nr) = AR COmer ICTENTs 
VAR - DRIVING NOISE VARIATION 
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CG?) C2 0) (2°) C2) GC) 


SUBROUTINE BURG(X,N,A,IP, VAR) 
COMPLEX X(N) ,ACO:N) ,EFK(500) ,EBK(500) 

COMPLEX EFK1(500) ,EBK1(500) ,AA( 20,20) ,SUMN,SUMD 
REAL RHO(0: 80) ,FPE(0: 80) 

INTEGER START 

RHO(0)=0. 

A(0)=CMPLX(1. ,0. ) 


DO 


DO 


DO 


10 I=1,N 
RHO( 0)=RHO( 0)+CABS(X(I) )**2/FLOAT(N) 
20 I=2,N 
EFK1(1I)=X(1) 
EBRK1( 1-1)=X(1-1) 
1000 K=1,IP 
SUMN=CMPLX(0. ,0. ) 
SUMD=CMPLX(0. ,0. ) 
DO 30 I=K+1,N 
SUMN=SUMN+EFK1(1I)**CONJG(EBK1(I-1)) 
SUMD=SUMD+CABS(EFK1( 1) )**2+CABS( EBK1( 1-1) )**2 
AACK ,K)=-2. *SUMN/SUMD 
TENP=1. -CABS(ASGIn 
IF(TEMP. LE. 0.) TEMP=1. E-10 
RHO( K)=TEMP**RHO( K-1) 
IF(K.GT.1) THEN 
DO 40 J=1,K-1 
AA(J,K)=AA(J,K-1)+AACK,K)*CONJG( AA(K-J,K-1)) 
ENDIF 
DO 60 I=K+2,N 
EFK( 1)=EFK1(1)+AA(K,K)*EBK1( 1-1) 
EBK( I-1)=EBK1(1-2)+CONJG( AA(K,K) )*EFK1(I-1) 
DO 70 I=K+2,N 
EFK1(1)=EFK(I) 
EBK1(I-1)=EBK(I-1) 


CONTINUE 
BO 100" 1=i. Le 


ACI)=AACI,IP) 


VAR=RHO( IP) 
RETURN 


END 


CHAN SEa EE. 


PEO T 

B - AR COEFFICIENTS 

A - TEMPORARY USING Saka 

ISIZE - ALOG(CFFT DATA POINTS) /ALOG2 
OUTEUL 

A - POWER SPECTRUM 


eae@epeweegpeeeeeweeqew»eweueaqeeeawvweeedqrFb@@3 @]ePerwerweweweswwFeswesw ws seg BPe*s28 Fw Bw BT FT BT FZ Te Fs TF se FW wZwF82 SF FBT FT FXx Fie #e#*#e#Oeog_ZLe28reeee2 8 FF BF Wi 
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CoC) Cy GPa a ae 


Cay Ce CoCo 0) C2 <7 


200 


SUBROUTINE CHANGEFFT(A,B,ISIZE,K) 
COMPLEX A(0: 2°"ISIZE-1) ,B(0: 2**ISIZE-1) 
INTEGER REVRSE 
PI=ACOS(-1. ) 
N=2*"ISIZE-1 
K=0 
FS=FLOAT(N+1) 
Chimeric N Size A. bons) 
CALL FINDMAX(A,N,K) 


RETURN 
END 
FINDMAX 
INPUT 
N - NUMBER OF SPECTRUM POINTS(N+1) 
A ="POWBR SPECTRUM 
OUTPUT 
MAX - ARRAY INDEX OF MAXIMUM POWER SPECTRUM 


eeeeeeeeweeeweeeewenrnewereewenereereeweneweeeeweweweweeweweweweweeeweweweeweweweweweeweweweweweeweweweeweweweweweewee@ @ @ @ ow %. 


SUBROUTINE FINDMAX(A,N,MAX) 

COMPLEX A(0:N) 

MAX=0 

AMAX=1. E66 

DO 100 I=0,N 
IF(CABS(A(1I)). LT. AMAX) THEN 


MAX=I 
AMAX=CABS( ACI) ) 
EN DEE 
CONTINUE 
RETURN 
END 
REVERSE 
BIT REVERSE ORDER INDEX CHANGING SUBROUTINE 
INPUT 
N - ARRAY INDEX 
ISIZE - ALOG(FFT DATA POINTS) /ALOG2 
OOTEUT 


REVERSE=) Bil REVERSE ORDER INDEX 


INTEGER FUNCTION REVRSE(N,ISIZE) 
INTEGER A(20) 
DO 200 I=1,ISIZE 
A(1)=MOD(N,2) 
N=N/2 
CONTINUE 
REVRSE=0 
DO 300 I=1,ISIZE 
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200 


300 


400 
500 


REVRSE=REVRSETA( 1 27 Ga ee 


CONTINUE 
RETURN 
END 
DFT 
COOLEY - TUKEY ALGORITHM 
INPUT 
N - NUMBER OF DATA POINTS 
ISIZE - ALOG(N)/ALOG2 
A - BIT REVERSE ORDER INDEXED TIME DOMAIN 
B - TEMPORARY ARRAY 
CUPUL 
A - FREQUENCY DOMAIN 


eee eweeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeww e@e ew e2eene ee 2 @ @ @ = %* 


SUBROUTINE DFT(N,ISIZE,A,B) 
COMPLEX ACO: je BCC 
PI=ACCSG- 13) 
DO 500 -t1= size 
ISAGE1=2***(I1-1) 
ISTAGE=2***T 1 
DO 200 I=0,N 
TTEST=NODC 1. PS TAGE) 
IFCITEST. GT. ISAGE1) THEN 
K=ITEST-ISAGE1 
THEATA=2. *PI*FLOAT(K)/FLOAT( ISTAGE ) 
T1=SINCTHEATA) 
T2=COS( THEATA) 
TECABS(UT 1) 2br 15 => elo 
LFCABS( 12)... 1E=5 ) gZ2—0 
W=CMPLX(T2,-T1) 
ACI)=AC1)*W 
ENDIF 
CONTINUE 
DO 300 I=0,N 
ITEST=JOD(2 , 1S TAGE) 
TFCITES?. GE. 1SAGEI > THEN 
BC I)=-AC I )+AC I-ISAGE1 ) 
Buod 
BCI)=AC1)+AC I+ISAGE1) 
ENE LF 
CONTINUE 
DO 400 I=0,N 
AC 1I)=BCI) 
CONTINUE 
CONTINUE 
RETURN 
END 
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CC) CG Gea a CaCl (a GG 2) 


1. INDEX CHANGE 
Z2.USE COOLEY TUKEY ALGORI TM 


INPUT 
N - NUMBER OF DATA POINTS 
ISIZE - ALOG(N)/ALOG2 
B - INPUT(TIME DOMAIN) 
FS - SAMPLING FREQUENCY 
OUTEUT 
A - FREQUENCY DOMAIN 


Sen ROUMINE FRICN,ISIZE,A,B,FS) 
COMPLEX A(0:N),BC0:N) 
INTEGER REVRSE 
DO 100 I=0,N 
k=I 


eee evwvnwrnwrwroeeeaeaeaeeeeaenweaeweeeaeeeeweaeaeeeweaenwaeewaeewewenwraeewaeewnwreeweweeweeweeweeweeeweweweeeaeeweee ws. 


100 ACREVRSE(K, ISIZE) )=BC 1) 
CAGL Diane SIZE, AsB) 
DO 200 I=0,N 
200 AC I)=AC1)/FS 
RETURN 
END 
C iain Rentart act asian tant aslastastart acfast ant ari astasiasiasiast axtaniastantastastastasiastastasiaslasiasiasiasiasiasiastasianiastas| an anjanas/as/es as as anes as! as an as as as an on am as os on oo 
e 
v CLEAR 
C NEC 
C N - NUMBER OF DATA POINTS 
C A,B - INPUT arrays 
C OULeUT 
C A,B - reinitiallized arrays 
C 
ee 
SUBROUTINE CLEAR(A,B,N) 
COMPLEX ACO: N),BC 0: N) 
DO 100 I=0,N 
100 AC I)=B( I)=CMPLX( 0. ,0. ) 


RETURN 
END 
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